function [Y,XX,W] = plain(S,IsDeviation,A0,Ea,Eu,NPer,V)
% simulate.linear.plain  [Not a public function] Plain linear simulation.
%
% Backed IRIS function.
% No help provided.

% -IRIS Toolbox.
% -Copyright (c) 2007-2015 IRIS Solutions Team.

%#ok<*VUNUS>
%#ok<*CTCH>

% The input struct S must at least include
%
% * First-order system matrices: T, R, K, Z, H, D
% * Effect of nonlinear equations: Q, v
%

try
    V;
catch
    V = [ ];
end

%--------------------------------------------------------------------------

% First-order solution matrices.
T = S.T;
R = S.R;
K = S.K;
Z = S.Z;
H = S.H;
D = S.D;

ny = size(Z,1);
nxx = size(T,1);
nb = size(T,2);
nf = nxx - nb;
ne = size(Ea,1);
R0 = R(:,1:ne);
colR = size(R,2);
isShkSparse = issparse(Ea);

if IsDeviation
    K(:) = 0;
    D(:) = 0;
end

Y = nan(ny,NPer);
W = nan(nxx,NPer); % W := [xf;alp].

lastEa = max([ 0, find(any(Ea ~= 0,1),1,'last') ]);
lastEu = max([ 0, find(any(Eu ~= 0,1),1,'last') ]);

% Nonlinear add-factors.
IsNonlin = ~isempty(V) && ~isempty(S.Q);
if IsNonlin
    Q = S.Q;
    lastN = utils.findlast(V);
    colQ = size(Q,2);
else
    lastN = 0;
end

% Initial condition.
if isempty(A0)
    wt = zeros(nxx,1);
else
    wt = [ zeros(nf,1); A0 ];
end

% Transition variables
%----------------------
for t = 1 : NPer
    wt = T*wt(nf+1:end) + K;
    if lastEa>=t
        eat = Ea(:,t:lastEa);
        eat = eat(:);
        nAdd = colR - size(eat,1);
        if isShkSparse
            eat = [ eat; sparse(nAdd,1) ]; %#ok<AGROW>
        else
            eat = [ eat; zeros(nAdd,1) ]; %#ok<AGROW>
        end
        wt = wt + R*eat;
    end
    if lastEu>=t
        wt = wt + R0*Eu(:,t);
    end
    if lastN>=t
        vt = V(:,t:lastN);
        vt = vt(:);
        nAdd = colQ - size(vt,1);
        vt = [ vt; zeros(nAdd,1) ]; %#ok<AGROW>
        wt = wt + Q*vt;
    end
    W(:,t) = wt;
end


% Mesurement variables
%----------------------
if ny>0
    Y = Z*W(nf+1:end,1:NPer);
    if ~isempty(Ea)
        Y(:,1:lastEa) = Y(:,1:lastEa) + H*Ea(:,1:lastEa);
    end
    if ~isempty(Eu)
        Y(:,1:lastEu) = Y(:,1:lastEu) + H*Eu(:,1:lastEu);
    end
    if ~IsDeviation && any(D(:) ~= 0)
        Y = Y + repmat(D,1,NPer);
    end
end

XX = [ W(1:nf,:); S.U*W(nf+1:end,:) ];

end
